%solving strategy using residual based incremental update static scheme
%with combination with residual based newton raphson strategy
function [elements,output] = solving_strategy_static_1(elements,model,...
    numTimeStep,maxits,tol,postProcessInfo)
% tol = 1e-6;
% maxits = 50;
dim = 2;
npp = size(postProcessInfo.infoList,2);
output = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time steps iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.preAnalysis();
mesh = model.getMesh();
nDofs = mesh.countDofs();
r = zeros(nDofs,1);
r = model.assembleR(r);
u = zeros(nDofs,1);
lambda = 0;
dlambda = 1/numTimeStep;
for i=1:numTimeStep
    lambda = lambda+dlambda;
    lr = lambda*r;
    %update model's states
    model.setX(lambda,u);
    elements = InitializeSolutionStep(elements,model);
    
    %assemble ri
    ri = assembleRi(elements,model);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Newton-Raphson loop
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cnt = 0;
    while 1
        cnt = cnt+1;
        disp(['substep ' num2str(cnt) '--------------']);
        %assemble kt
        kt = assembleKt(elements,model);
        
        du = kt\(lr-ri);
        u = u + du;
        model.setX(lambda,u);
        
        %update model's states
        elements = CalculateMaterialResponse(elements,model);
        
        %assemble ri
        ri = assembleRi(elements,model);
        
        %compute convergent criteria
        %res = norm(du)/norm(u);
        res = norm(r-ri);
        
        disp(['substep ' num2str(cnt) ' complete']);
        if res < tol
            disp(['Newton Raphson iteration converged in ' ...
                num2str(cnt) ' loops']);
            break;
        end
        if cnt>maxits
            disp(['Newton Raphson iteration did not converge in ' ...
                num2str(maxits) ' loops']);
            break;
        end
    end
    
    elements = FinalizeSolutionStep(elements,model);
    
    %     findmodelaximumPrincipalStress(elements);
    %
    %     responsePlot.x(i) = i*ddu;
    %     responsePlot.y(i) = elements(398).stress(1);

    output = processOutput(i,dim,elements,model,...
        output,postProcessInfo);
    
%     for j=1:npp
%         req = postProcessInfo.infoList{j};
%         if strcmp(req,'resistant force')
%             riFull = computeResistantForce(elements,model);
%             param = postProcessInfo.infoParam{j};
%             nn = size(param,2);
%             cnt2 = 1;
%             for k=1:nn
%                 nodeid = param(k);
%                 output.resistant_force(:,cnt2,i) = ...
%                     riFull([nodeid*dim-1 nodeid*dim]);
%                 cnt2 = cnt2+1;
%             end
%         end
%         if strcmp(req,'stress at Gauss points')
%             param = postProcessInfo.infoParam{j};
%             ne = size(param,2);
%             for k=1:ne
%                 eid = param(k);
%                 output.stress_gp(:,:,k,i) = ...
%                     elements(eid).stress;
%             end
%         end
%     end
    
    disp(['time step ' num2str(i) ...
        ' finished with lambda = ' num2str(lambda)]);
end

end
